%%%%%%%%%%%%%%%%%%
% This is a TEMPLATE script for running the OHBA recommended 
% pipeline for doing source space ERF OAT analysis
% on Elekta-Neuromag data (a very similar pipeline will work on
% CTF data as well). 
%
% This takes as input data that has been preprocessed (e.g. using the OHBA
% recommended prepocessing pipeline (see the OSL wiki).


%%%%%%%%%%%%%%%%%%
%% SETUP THE MATLAB PATHS
% make sure that fieldtrip and spm are not in your matlab path

global OSLDIR;
    
%tilde='/home/mwoolrich/Desktop';
tilde='/Users/woolrich';
osldir=[tilde '/homedir/matlab/osl1.3.1'];    

addpath(osldir);
osl_startup(osldir);

%%%%%%%%%%%%%%%%%%
%% SPECIFY DIRS FOR THIS ANALYSIS

% directory where the data is:
datadir=[tilde '/homedir/vols_data/murphy_project/data']; % this is the directory the oat files will be stored in

% this is the directory the preprocessed SPM files are stored in:
spmfilesdir=[datadir '/spm_files']; 

% this is the directory the analysis files will be stored in:
workingdir=[datadir '/results']; 
cmd = ['mkdir ' workingdir]; unix(cmd); % make dir to put the results in

%%%%%%%%%%%%%%%%%%
%% Set up the list of subjects and their structural scans for the analysis 

clear spm_files structural_files

% Setup a list of SPM MEEG object file names to be created, in the same order as spm_files and fif_files:
% Note that here we only have 1 subject, but more generally there would be
% more than one, e.g.:
% spm_files{1}=[workingdir '/spm8_meg1.mat'];
% spm_files{2}=[workingdir '/spm8_meg1.mat'];
% etc...
is=15:28;
for i=1:length(is),
    spm_files{i}=[spmfilesdir '/efspm8_nosss_meg' num2str(i) '_1.mat'];
end;

% Setup a list of existing structural files, in the same order as spm_files and fif_files:
% e.g.:
% structural_files{1}=[datadir '/structurals/struct1.nii'];
% structural_files{2}=[datadir '/structurals/struct1.nii'];
% etc...
structuralsdir=[datadir '/structurals']; % this is the directory the structural files will be found in
% structural files:
sfs={'F9'    'F11'    'F15'    'F16'    'F20'    'F22'    'F25'    'M3'    'M6'    'M7'    'M12'    'M14'    'M18'    'M19'    'F4'    'F7'    'F10'    'F13'    'F18'    'F24'    'F26'    'M2'    'M4'    'M8'    'M10'    'M13'   'M17'    'M21'    'F3'    'F6'    'F12'    'F14'    'F17'    'F19'    'F23'    'M5'    'M9'    'M11'    'M15'    'M20'    'M22'    'M23'};
sfs{21}=''; % structural is missing
sfs{14}=''; % structural is missing
sfs{6}=''; % structural is missing
for i=1:length(is),
    if ~strcmp(sfs{is(i)},'')
        structural_files{i}=[structuralsdir '/' sfs{is(i)} '/struct.nii'];
    else
        structural_files{i}='';
    end;
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% setup oat

oat=[];

%%%%%%%
% source recon settings
oat.source_recon.D_continuous=[]; % do not have continuous files, but if you do then it is best to specify these too
oat.source_recon.D_epoched=spm_files; % epoched files
oat.source_recon.conditions={'Motorbike','Neutral face','Happy face','Fearful face'};
oat.source_recon.freq_range=[1 48]; % frequency range in Hz
oat.source_recon.time_range=[-0.2 0.5];
oat.source_recon.method='beamform'; % switch this to 'none' to do sensor space analysis
oat.source_recon.gridstep=8; % in mm, should be >=6mm
oat.source_recon.mri=structural_files;
oat.source_recon.sessions_to_do=1:14; % sessions (in this case subjects) to run over
oat.source_recon.dirname=[workingdir]; % where the oat directory will be saved

%%%%%%%
% first level settings

% setup design matrix (see OSL wiki for more on this)
Xsummary={};Xsummary{1}=[1 0 0 0];Xsummary{2}=[0 1 0 0];Xsummary{3}=[0 0 1 0];Xsummary{4}=[0 0 0 1];
oat.first_level.design_matrix_summary=Xsummary;

% contrasts to be calculated:
oat.first_level.contrast={};
oat.first_level.contrast{1}=[3 0 0 0]'; % motorbikes
oat.first_level.contrast{2}=[0 1 1 1]'; % faces
oat.first_level.contrast{3}=[-3 1 1 1]'; % faces-motorbikes
oat.first_level.contrast{4}=[0 0 -1 1]'; 
oat.first_level.contrast{5}=[0 -1 0 1]';
oat.first_level.contrast{6}=[3 1 1 1]';
oat.first_level.contrast_name{1}='motorbikes';
oat.first_level.contrast_name{2}='faces';
oat.first_level.contrast_name{3}='faces-motorbikes';
oat.first_level.contrast_name{4}='fear-happy';
oat.first_level.contrast_name{5}='fear-neutral';
oat.first_level.contrast_name{6}='faces+motorbikes';

%%%%%%%
% group-level stage
%
% No settings specified at all - so this will just use defaults to do a simple group
% average.
% See OSL wiki for how to setup group design matrices for other scenarios

%%%%%%%
% call func to set the settings
oat = osl_check_oat(oat);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Run the oat

oat.to_do=[1 1 1 1]; % run ALL stages 

oat = osl_run_oat(oat);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% OUTPUT INDIVIDUAL SUBJECT NIFTII FILES

str=['fslview ' OSLDIR '/std_masks/MNI152_T1_' num2str(oat.source_recon.gridstep) 'mm_brain '];
for ii=1:1,%14,
    S2=[];
    S2.oat=oat;
    S2.stats_fname=oat.first_level.results_fnames{ii};
    S2.first_level_contrasts=[3]; % list of first level contrasts to output
    S2.resamp_gridstep=oat.source_recon.gridstep;

    [statsdir,times]=osl_save_nii_stats(S2);
    str=[str statsdir '/tstat3_10mm '];
end;

str2=[str ' &'];
runcmd(str2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% OUTPUT GROUP NIFTII FILES

S2=[];
S2.oat=oat;
S2.stats_fname=oat.group_level.results_fnames;
S2.first_level_contrasts=[1:3]; % list of first level contrasts to output

[statsdir,times]=osl_save_nii_stats(S2);
